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ABSTRACT 

Millisecond pulsars are believed to be old pulsars spun up by a surrounding accretion disc. Magnetic fields are thought to play a 
leading role in this, both by determining the location of the inner edge of the disc and by exerting an additional torque on the star (as 
a result of the interaction between the stellar magnetic field and the disc plasma motion, which creates a toroidal component B^). In 
some well-known analytic models, developed in the 1980s, the profile was taken to be proportional to the relative angular velocity 
between the disc plasma and the neutron star, multiplied by a vertical dipolar field. The present work stands in the line of improving 
those models, suggesting a new profile for B. In a previous paper, we discussed the poloidal component of the magnetic field and 
here we consider the toroidal component, again making the kinematic approximation and looking for steady solutions of the induction 
equation for axisymmetric models. The poloidal magnetic field is not assumed to be dipolar and the poloidal velocity field is not 
taken to be zero everywhere. We also do not use the thin disc approximation to simplify the induction equation but instead solve it 
numerically in full 2D. The profile obtained in the earlier analytic models is shown to have very limited validity and a more general 
semi-analytic solution is proposed. 

Key words, accretion, accretion disks - magnetic fields - magnetohydrodynamics (MHD) - turbulence - methods: numerical - 
X-rays: binaries 



1. Introduction 

In the present work we study the deformation caused in a neu- 

■ tron star's magnetic field because of the interaction with the mat- 

■ ter in a surrounding accretion disc. A basic description for this 
kind of system was given by Ghosh and Lamb in 1979j with 
the model subsequently being improved by Wang (1987]l and 

' Campbell (I1987I I. who suggested an analytic expression for the 
, toroidal component of the field. This expression has been widely 
used since then, on account of it being both simple and physi- 
cally plausible. 

In these analytic models, the authors made the kinematic ap- 
' proximation, looking for an axisymmetric stationary solution of 
' the induction equation with a given unchanging structure for 
the fluid in the disc. They further took the disc to be thin, the 
poloidal component of the magnetic field to be exactly dipolar, 
and the velocity field to have zero poloidal component, with its 
azimuthal component being Keplerian inside the discQ and coro- 
tating at the top of the disc. Using cylindrical coordinates (m, (f>, 
z), they found 

B^^y^(nK-n.)B,TiCcAQ./m\ (1) 

where is the amplification factor, D.k and are the Keplerian 
and stellar angular velocities, respectively, and is the dissipa- 

' Campbell also considered non-Keplerian flow in the inner part of 
the disc. 



tion time scale. The amplification factor was taken to be a 
constant not much greater than 1 (it depends on the steepness 
of the transition - in the z-direction - between Keplerian mo- 
tion inside the disc and corotation with the star at the top of the 
disc). The precise profile of depends on what is the domi- 
nant mechanism for dissipating the magnetic field. Wang (119951) 
considered three different cases, with Td being dominated by the 
Alfven velocity, turbulent diffusion and magnetic reconnection, 
respectively. 

Equation ^ was then used for calculating the net magnetic 
torque exerted on the neutron star. The vertically averaged torque 
can be written as 

Tec B^B- m/h , (2) 

where h is the semi-thickness of the disc. Regions of the disc in- 
ward of the corotation point therefore give positive contributions 
to the torque (because B^ > 0), while the remainder of the disc 
gives negative contributions (because B^ < 0). The total mag- 
netic torque is obtained by integrating the local values from the 
inner edge to the outer boundary, and it can be either positive or 
negative depending on the location of the inner edge of the disc 
with respect to the corotation point. 

The aim of the present paper is to develop a semi-analytic 
model that can improve on those of Wang and Campbell, while 
remaining simple enough to be useful for people discussing the 
behaviour of astrophysical sources, giving a conceptual picture 
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to go alongside results from large numerical calculations where 
the full set of the MHD equations is solved. 




ourer aisc 
(main region) 



Fig. 1: Schematic representation of our model (not to scale). We 
use Tin = 10 Tg, rti- ~ 22 and Hc ~ 115 r^. The opening angles 
are 8° for the disc alone and 10° for the disc plus corona. The 
outer disc extends much further out than the main region shown 
here: the grid continues until ro^t - 380 rg. 



Our approach, for the time being, is to continue to retain 
axisymmetry and the kinematic approximation but to calculate 
a consistent steady-state solution for the magnetic field, relax- 
ing the assumptions on the poloidal components of the magnetic 
and velocity fields and using a 2D model without any vertical 
averaging of the Taylor expansion of the induction equation. In 
the main region of the outer disc (see Fig. [TJ we use a simple 
Keplerian velocity profile, but this is something that will be im- 
proved on later In a previous paper (Naso & Miller l 20T0l here- 
after Paper I) we analysed the distortion of the poloidal compo- 
nent of the magnetic field using a similar approach, and found 
that deviations away from a dipole field can be quite significant. 
Here we focus on the toroidal component and use the results of 
the previous model to solve the cp component of the induction 
equation. We find that in general follows a profile different 
from that of the analytic models, i.e. Eq. ([T]|, and reduces to that 
only in a very particular case. 

Following this introduction, in Sect.|2]we briefly describe our 
model, which is the same as that of Paper I; in Sect. [3] we recall 
the equations used (obtained from the induction equation), give 
expressions for the velocity and diffusivity profiles and outline 
our solution method (details of tests made on the code are given 
in an Appendix); in Sect. |4] we present our numerical results; in 
Sect. |5] we comment on these, comparing them with those of the 
earlier analytic models, and develop our new suggestion for the 

profile. Sect.|6]contains conclusions. 



2. Model 

In this study, we are considering disc accretion by a neutron star 
having a dipolar magnetic field. The model is the same as that 
of Paper I. For a detailed description of it, see Section 2 of that 
paper; here we recall the main points. 

We are assuming that the star is rotating about its magnetic 
axis, and that this axis is perpendicular to the plane of the disc; 
also, we assume that the fluid flow is steady and has axial sym- 
metry everywhere. We use the kinematic approximation and do 
not consider any dynamo action, but turbulent diffusivity is in- 
cluded. The velocity field is not constrained to be purely az- 
imuthal but is allowed to have non-zero components also in the 
other directions. We use spherical coordinates (r, 0, (p), with the 
origin being at the centre of the neutron star. Boundary condi- 
tions are imposed at the inner and outer radial edges of the disc 
(rin is at the Alfven radius, and rout is at 38 ^n), on top of the 
corona (taken as being a layer above and below the disc) and 



on the equatorial plane. Having the inner edge of the disc at the 
Alfven radius justifies the kinematic approximation to some ex- 
tent, since in this configuration the magnetic pressure is smaller 
than the gas pressure within the region that we are considering, 
and so the effects of the plasma on the magnetic field should be 
larger than the magnetic feedback on the plasma flow. The ratio 
h/r is taken to be constant, with the opening angle being 8° for 
the disc (measuring from the equatorial plane to the top of the 
disc), and 10° for the disc plus corona. 

The inner disc region (r < r^: ~ 2 ri,,) and the corona are 
modelled with a larger value of 77 than the other parts. In these 
regions, the kinematic approximation does not provide a good 
description of the system for two different reasons: in the corona 
this is because of the low density of the plasma (which there- 
fore tends to follow the magnetic field behaviour rather than be- 
ing followed by it); in the inner region, it is because the mag- 
netic field intensity is still quite large - although the magnetic 
pressure is smaller than the gas pressure, it is not yet negligible. 
Using a larger value for in these regions makes the magnetic 
field less sensitive to the plasma motion; a somewhat similar 
approach was used by Kueker et al. (2000). We recall that the 
present knowledge of the turbulent magnetic diffusivity is quite 
poor and it is not a simple task to obtain a reliable expression for 
the 77 profile. 

As regards the velocity field: for v,. we use the expression 
given for the "middle region" of a-discs by Shakura & Sunyaev 
(I1973I I. For Q. we take Keplerian rotation in the disc and corota- 
tion at the top of the corona and at the inner edge of the disc, giv- 
ing a maximum for Q. between and r^j. These different parts 
are smoothly connected using error functions. Regarding vg: we 
put it to zero in the disc but near to the boundaries we are forced 
to have a non-zero value in order to be consistent with the dipo- 
lar boundary conditions (as shown in Section 3.2 of Paper I) and 
so we use a non-zero profile in the corona. In this way we are in- 
cluding in the model an outflow from the surface of the corona, 
and this is in agreement with recent hydrodynamic simulations 
of accretion flows (Jiao and Wu. l20I fb . 

Summarising, we divide the surroundings of the central ob- 
ject into four parts (see Fig. [T] which is repeated from Paper I): 
(1) the inner disc, extending from ri,, out to a transition radius 
rtr ~ 2 rin (where the diffusivity changes); (2) the outer disc, ex- 
tending from rtr to an outer radius rout = 38 rin; (3) a corona, 
above and below the disc; and (4) everything else, which we 
take here to be vacuum. As a unit for radial distances, we use 
the Schwarzschild radius rg. Within the outer disc, we focus on 
what we call the main region, extending from rtr out to the light 
cylinder at ric ~ 1 1 rin. 



3. Equations 

In the kinematic approximation, one assumes that the velocity 
field remains fixed as specified, and the interaction between the 
magnetic field and the plasma is then described by the induction 
equation alone. In the presence of turbulence, it is more conve- 
nient to write this equation for mean fields rather than for the 
actual fields (which contain fluctuating parts as well). 
The time dependence of the mean field is given by 



5,B = V X (v X B + £ - 77ohmV x B) 



(3) 



where T/ohm = c^/4-ncr is the Ohmic diffusivity and S is the tur- 
bulent electromotive force. A common procedure is to expand £ 
in terms of the mean field and its derivatives and within the first 
order smoothing approximation one has fi - ajB - rijV x B, 
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where the ajB term generates the so-called a-effect. As in Paper 
I, we are neglecting this effect here and the induction equation 
then reduces to 



«9,B = V X (v X B - //V X B) , 



(4) 



where rj - //ohm + ™d is ~ t/t, because the turbulent difFu- 
sivity is much stronger than the Ohmic one. 

We note that the effects of a dynamo action on the disc struc- 
ture have recently been studied by Tessema & Torkelsson (I2010al 
and l2010 F). who estimated the toroidal magnetic field generated 
by the dynamo to be about an order of magnitude larger than 
the calculated according to the early models. Here we show 
that the profile of the toroidal field can be very different from the 
one suggested by those models, if the poloidal component is not 
forced to be a dipole, but is instead calculated self-consistently. 

As described in Paper I, our strategy consists of writing Eq. 
(|4| in spherical coordinates, applying the axisymmetry and sta- 
tionarity assumptions (i.e. putting d^l. . .] - . .] = 0) and then 
solving the final equations with the velocity field and magnetic 
turbulence profile given by the model. The three components of 
the induction equation are 



= dg lysine 
= 5r|r 
^ drir 



v,.Be - veBr - - [dr(rBe) - deB,-] 
r 

VrBg - VgBr - -[dr(rBg) - dgBr] 



r |r V^Br - VrB^ + ^dr(rB^) 
dg I vgB^ - v^Bg l—:dg{B^ sin 0) 



(5) 
(6) 

(7) 



The first two equations contain only poloidal quantities and have 
been solved in Paper I (making use of the magnetic stream func- 
tion). Here we focus on the third equation and solve it using the 
results for B^ and Bg from the previous analysis. 

We rewrite Eq. ^ in the following dimensionless way: 



dlBs + agg dlB0 + a^ d_,B^ + ag dgB^ + ai -i- ao = , 



(8) 



where x = r/rg (not to be confused with the Cartesian coordinate 
X used for the plots), B^ is the toroidal field in units of a reference 
field (as described in Sect. 13.41 below) and the dimensionless a 
coefficients are 



ag = 



1 

dgJ] 
x^rj 



1 COS0 



XT] 



dgTj 1 COS 6 

XT] X sinO 



1 



ai 



d,T] ^ 2 Vrrg 

T] X T] 

d,,T] rgdx(xvr) dgTj cos 
XT] XT] x^r] sin0 

1 cos^e r„dgvg 1 



(9) 
(10) 

(11) 



sin^ 



1 



XT] X'^ 
rgdxixVr) dgTj COS 

XT] sm0 



(12) 



1 COS^ r„dgVg 



X sin^ 



1 



ao 



rgdjcixv^Br) rgdgiv^Be) 



XT] 



XT] 



(13) 



All of these coefficients have direct analytic expressions, ex- 
cept for the last one (ao) which contains B, and Bg, whose values 
are taken from the previous numerical calculations. Note that the 
magnetic field components appearing in Eq. ( fT3T l are dimension- 
less. 



3.1. Poloidal velocity and diffusivity 

The poloidal components of the velocity field and the diffusivity 
have already been discussed in Paper I. Here we use the same 
profiles as before; for the sake of clarity, we give the expressions 
again below. 
For Vr, we use 



2 X 10^ a*^^in^'^m 



-1/5 



x{3/xf"[l-(3/xf'Y^" [cms-i], 



(14) 



where m is the stellar mass in solar mass units, m is the accre- 
tion flux in units of the critical Eddington rate and a is the stan- 
dard Shakura-Sunyaev viscosity coefficient. Using typical values 
(a - Q.l, m - 0.03 and m - 1.4) this gives 



vv(x) * 7.3 10^ ■ (3/xf'^ [l - (3/xf^] ^'^ cms-' 
For vg, we use 



ve(r,0) 



in the disc 

I 



(15) 



(16) 



^ V, tan in the corona 
where the transition in vg between the two regions is made using 



ve(r,0)^fi(0)^v, tan0 



with 



1 + erf 



(^) 



(17) 



(18) 



where 0b is at the upper surface of the disc (i.e. 82°) and d 
10-^ radians (i.e. 0.057°). 
For the diffusivity, we use 



T](r,0) = T]o\l+[T]g(0) + T]Ar)] 



^-1 



(19) 



where r]o and tj^ are the values in the main disc region and the 
corona (see Fig. [T]), for which we use the values lO'" cm^ s ' 
and lO'^ cm^ s"' respectively (we also use other values to study 
the impact on the results of varying 77). For T]g{0) and T]rir) we 
use joining functions of the form 



fix) = - 



1 -I- erf 



/ —X + Xc \ 

\ d / 



(20) 



with Xc - 0D and rt, for rjg and 77,., respectively, and with d being 
the width of the transition in the error function, for which we use 
J = 5 r„ in the radial direction and li = 2° in the directior0. 



^ Note that the error-function widths in the radial and angular direc- 
tions are larger than those used in Paper I (there = 2r^ and dg = 10"^^ 
radians = 0.057°). The reasons for this are explained in Appendix I A.2I 
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3.2. Azimuthal velocity 

We are modelling the main disc region as rotating with Keplerian 
velocity, with the corona being taken as a transition layer where 
the velocity goes from Keplerian to corotation. In the inner re- 
gion, we join the Keplerian flow to corotation at the radial inner 
edge, again using an error function. We recall that in Paper I we 
showed that a strictly dipolar field, without distortions, can in 
principle be a stationary solution of the induction equation pro- 
vided that the velocity field fulfils the two conditions: 



1 

Vg = — tanOvr 
and 

Q cc r^'^ sin'' . 



(21) 
(22) 



From Eq. (l22l l one sees that corotation, which is obtained by 
choosing y = 0, is consistent with having dipolar conditions 
(which is what we are using here). 

For the magnetic field intensities and neutron star spin rates 
which we are using as standards {B ~ 10^ G and P ~ 10 ms), 
the corotation point is outward of the inner edge of the disc 
(which is the standard condition required for being in the ac- 
cretion regime) and therefore Q. should reach a maximum at 
some location and then decrease again, as one moves inwards. 
Summarising, we use the following profile: 



Q(r, 0) 



smooth join in 
smooth join in r 



in the main region: 

0e [eD,7T/2],re[rn-,r,, 

in the corona: 

0e[ec,0D] 

in the inner disc: 
r e [rin, r^] 

at ghost zones: - 0c - A0 
at inner edge: r - 



(23) 



where 0c is the upper surface of the corona, /S.0 is the angular 
resolution, Qs is the stellar spin rate, D.k is the Keplerian angu- 
lar velocity ^jGM|^' and the two smooth connections are made 
using the error functions given in Eqs. (|25] | and (l27l i below. 
As regards the smooth joins, in the direction we write 



Q(r, 0) 



^K(r)fx(0) + ^Al-fxm 
[QHr)-Qs]/i(^?) + a 



where: 



Me) 



1 -Herf 



(^) 



(24) 



(25) 



with d - \Q ^ radians (i.e. 0.57°), and we then have a modifica- 
tion in the radial direction giving 



Q(r, 0) 



Q.{r.0)f2{r) + ^Al-f2{r)] 
[Q(r,0)-Qs]/2(r) + a, 



where 



fiir) 



1 + erf 



(26) 



(27) 



with d = Ivg. See Fig. |2]for a contour plot of Q(r, 0) (made in 
terms of Cartesian coordinates x and z.) 




Fig. 2: Contour plot of Q(r, 0). This is the same for all of the 
configurations. The straight white line indicates the boundary 
between the disc and the corona. The corotation point is at r = 
18.8 r„. 



3.3. Boundary conditions forB^ 

In our model we take the region outside of the disc and its asso- 
ciated corona to be vacuum, and suppose that there is no toroidal 
component of the magnetic field there (i.e. B^ - 0). We impose 
the same condition also at the equatorial plane because B^ has to 
be antisymmetric across it. 

When B^ = 0, Eq. (|7]i reduces to ao = 0, i.e. 



dAxv^Br) + deiv^Bg) = 



(28) 



at the top of the corona, at the inner edge of the disc and on the 
equatorial plane. We need our choice of the boundary conditions 
(i.e. purely dipolar field and corotation) to represent a solution of 
this equation. In Paper I, we showed that corotation is consistent 
with a pure dipolar field, as already mentioned in the previous 
subsection (see Eq. (l22ll). In addition we note here that, on the 
equatorial plane, Eq. (l28b is trivially satisfied because B^ - Q 
and both and Bg are symmetric with respect to this plane. 

3.4. Solution method 

In order to get to Eq. ([8]l we first expanded out all of the deriva- 
tives present in Eq. O and then put the result into a dimension- 
less form using x = r/r„ as the radial coordinate and measuring 
B^ in units of B^, which is a reference value for the magnetic 

field, taken to be Bq, where Bo is the magnitude of the dipolar 
field at the stellar equator and ro is the dimensionless neutron- 
star radius. As in Paper I, we use canonical values for the mass 
and radius of the neutron star, 1.4 M© and 10 km respectively, 
and take Bq = 3 x 10^ G, as typical for a millisecond pulsar 

Equation ([S) is a non-homogeneous elliptic partial diff'eren- 
tial equation for B^ and we have solved it using the Gauss-Seidel 
relaxation method which uses a finite-difference technique, ap- 
proximating the operators by discretizing the functions over a 
grid. At any given iteration step, the values of B^ at the various 
grid points are written in terms of values at the previous step, or 
at the present step in the case of locations where it has akeady 
been updated. Details of the numerical scheme are given in the 
Appendix of Paper I, where we used the same method to solve 
the elliptic partial dififerential equation for the magnetic stream 
function. 

Before applying the numerical scheme to the actual problem 
that we want to solve, we performed a series of tests on the code. 
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which are described in detail in the Appendix. We used several 
configurations, with many different numbers of grid points, pro- 
files of the turbulent diffusivity, initial estimates for the toroidal 
field, locations for the outer radial boundary of the grid and val- 
ues for the iteration time step. In this way we have checked the 
code stability and convergence, have optimised the iteration pro- 
cedure and have determined where best to place the outer radial 
boundary of the grid (so that the outer boundary conditions do 
not significantly influence the solution in our region of interest). 

All of the results presented in this paper have been obtained 
using a grid spacing of Ar = 0.74 and A0 - 0.125°. The 
final maximum residual was of the order of 10"'^ and satura- 
tion was reached after about 4x10^ iterations, the iteration time 
step being 4.46 x lO"'* (about 95% of the critical one, beyond 
which the method gives a divergent solution). As a first estimate 
for the profile of the magnetic field we have used a Gaussian. 
The numerical domain was x € [10,380], 6 € [80,90] (in de- 
grees), which we covered with a homogeneous grid of 501 x 81 
points. The profile of the poloidal field (which is present in the 
expression for the coefficient oq) was calculated by running the 
code used in Paper I with the same profile of 77 and v and the 
same resolutions as used in this analysis. However, since for the 
poloidal calculation there is a stronger dependence on the outer 
radial boundary conditions, we have placed rout at 750 rg and 
used a 1001 x 81 grid. In general the number of radial points 
used in the poloidal calculation (A^f°') and in the toroidal one 
(Nf"^) are related through the following condition (which comes 
from equating the spatial resolutions): 



Nr - 1 = 



pol 
^out ~ 



(Nf-l). 



(29) 



For the poloidal calculation the final maximum residual was of 
the order of 10"^ and was reached after about 7 x 10'' iterations. 
The iteration time step was the same as for the toroidal calcula- 
tion. 

4. Results 

4.1. Reference configuration 

As mentioned in Sect. [3] we have used a slightly different con- 
figuration from that considered in Paper I (the transition in the 77 
profile is wider and the resolution in the angular direction is in- 
creased). The poloidal field for the new configuration, resulting 
from solving Eqs. (|5]l and (|6]l, is shown in Fig. [3] If we addition- 
ally assume a profile for the angular velocity, we can then solve 
Eq. (|8]l. The result, for the profile given by Eq. ( l23l l. is shown 
in Fig. |4]as a contour plot of B^, for all of the region of interest 
(i.e. for re [10, 1 15] r^ and including the corona). Contours for 
positive values of B^, representing toroidal field lines rotating 
in the same direction as the neutron star, are shown with black 
solid lines; while those for negative values of B^ (rotating in 
the opposite direction) are shown with black dashed lines. Triple 
dotted-dashed white lines show where B^ - 0. 

The toroidal field shows a quite structured profile, with two 
(positive) maxima and two (negative) minima. Surprisingly, the 
global maximum is located outward of the corotation point and 
is positive, in contrast to the standard picture according to which 
the sign of B^ is the same as that of the relative angular velocity 
between the disc matter and the star There is also a quite striking 
vertical structure, the global minimum being located just above 
and to the left of the global maximum. As in the analytic models, 
B^ becomes zero near to the corotation point (rcor - 18.8 rg), but 




Fig. 3: Poloidal field lines from the numerical solution (solid) 
and those for a dipole (dotted). In this configuration, vq - 10^ 



cm s while the diffusivity 770 



10^° cm2 s-i 



with a transition 



width in 9 direction of 2° (3.5 x 10"^ radians). Straight lines 
indicate the top surface of the corona (solid) and the boundary 
between disc and corona (dashed). 



73 0.91 




Fig. 4: Contour plot of the toroidal field in Cartesian coordinates, 
in the region r = [10 - 115] rg and = 80° - 90°. The black 
solid lines are lines of positive B^, while black dashed lines are 
used for negative B^. B^ - Q is shown with triple dotted-dashed 
white lines. The straight white line indicates the boundary be- 
tween disc and corona. 



here one sees that for almost every value of r > rcor there is a 




Fig. 5: The same as in Fig. |4]but showing an expanded view of 
theregionr = [10- 35] r^. The corotation point is at r = 18.8rg. 
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value of 6 where the field is zero, so that there are places with 
zero throughout all of the main disc region. 

In Fig. |5] we show the contour plot in the region r € 
[10, 35] Tg, so that the structure of the magnetic field in this part 
can be seen more clearly. In Table [T] we give the coordinates of 
the minima and maxima and the magnitude of the toroidal field 
at those locations, measured in units of the stellar field strength 
Bq (which is taken to be 3 x 10*^ G). All of these main features 
are located in the region r € [10,55] r^; the remainder of the 
disc can be divided into two zones: one where B^ > (in Figs. 
|4] and |5] this is below the long white triple dotted-dashed curve 
that crosses the equatorial plane at x = 23 and x = 1 10), and the 
other where, instead, it is negative. The minimum latitude of the 
region with positive B^ is about 84°. 



Table 1 : Locations of extrema of the toroidal field. 



Extremum 




e 


60 [5o] 


Global maximum 


30.5 


87.1° 


0.78 


Local maximum 


14.5 


83.1° 


0.092 


Global minimum 


25 


83.4° 


-0.17 


Local minimum 


52 


83.6° 


-0.13 



Notes. The last column reports the intensity of the toroidal magnetic 
field at the corresponding locations, measured in units of the stellar field 
strength Bq (which is taken to be 3 x 10^ G). 

Radial profiles of B^ are shown in Fig.|6]for several values 
of the latitude and for the shell average. We show the profiles 
for 6 = 87.1° (where the global maximum is), 6 = 83.6° (where 
the local minimum is, and which is very close also to the lo- 
cal maximum and the global minimum) and for an intermediate 
value - 85°. The curves all pass through zero at the corotation 
point (r - 18.8 r^), at least to within the accuracy of the calcu- 
lation. The large positive B^ peak at about 30 rg is progressively 
reduced as one moves from the mid plane of the disc towards the 
corona, and eventually becomes a negative local maximum. One 
can calculate the shell average of the radial profile, i.e. an aver- 
age of Z?0 over 6 for each r, and this is also shown in Fig. |6]as a 
thick solid line. This average reproduces quite well the general 
behaviour of the toroidal field and shows all of its main features. 
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Fig. 6: Toroidal field strength plotted against r at different val- 
ues of 0, in the region r - [10 - 50] rg (it is equal to zero on 
the equatorial plane). Negative values mean that B^ is pointing 
backwards with respect to the disc rotation. The thick solid line 
is the shell average of B^. 



The dependence of the sign of B^ is a key result, that can 
have quite dramatic consequences for the calculation of the mag- 
netic torque exerted on the neutron star In fact, it is usually as- 
sumed that the torque depends only on r, being positive inside 
the corotation radius and negative outside it (see Eqs. ([TJ and 
(IS)), whereas we now see regions of positive B^ even outward of 
the corotation point. Therefore we need to rethink the discussion 
of which regions in the disc tend to spin the star up or down. 
An appropriate approach for calculating the torque requires in- 
tegration in both directions: we plan to perform this analysis in 
a future work. 

Finally, we note here that the magnitude of the toroidal com- 
ponent of the magnetic field is typically larger than that of the 
poloidal component. For example, the maximum of the shell av- 
erage of the poloidal component + B] ^ jis ~ 3.6xlO-\ 

while the maximum of (^0)^ is ~ 2.3 x 10"'. In terms of en- 
ergy conservation, one has to bear in mind that any change in 
the magnetic energy has to be compensated by a corresponding 
opposite change in the plasma energy, while in the present model 
we are taking the flow pattern to be fixed. It is important for the 
back-reaction to be consistently taken into account and this will 
be done in subsequent work. Also, the distortions of the toroidal 
field will be limited by magnetic reconnection. 

4.2. Configurations witli largen]o 

We have already mentioned in Sect. |2]that there are big uncer- 
tainties about how to model the turbulent magnetic diffusivity 
within the disc and the surrounding corona. This quantity is of- 
ten discussed in terms of the turbulent magnetic Prandtl number, 

= v/r], which links it with the turbulent viscosity. 

Within the kinematic approximation, one does not solve self- 
consistently for the velocity field and so it is necessary assume 
some profile for it. As outlined in Sect. 13.1! in the present cal- 
culations we are using the velocity prescription given by the 
Shakura & Sunyaev thin disc model, which also embodies a par- 
ticular connection between the viscosity and other disc quan- 
tities. This gives (following e.g. Szuszkiewicz & Miller 120011) 
V = ff/iCs, with ^Ih^ = (assuming vertical hydrostatic 

equilibrium), where Cj is the sound speed, p is the pressure and 
p is the density. Using the isothermal sound speed, the Keplerian 
angular velocity profile and the parameters given in Sect. 13.1! 
one gets v ~ lO'-' cm^ s"' which, in turn, gives ~ lO-' in 
the disc and ~ 10 in the corona of our reference model. These 
values may seem rather high, but one should be cautious about 
using them because there are several further factors which need 
to be taken into account. 

Note, first, that calculating the viscosity with the a disc 
model certainly gives an over-estimate for P^, because one is ne- 
glecting the contribution of the magnetic field in the equation of 
motion. Also, there is a degeneracy in the profile of the velocity 
field, in the sense that one can obtain the same velocity profile, 
and hence the same results for our numerical calculations, using 
different combinations of a and m: an increase/decrease by a fac- 
tor of A in the accretion rate, together with a decrease/increase 
by a factor of yfJ in a (and hence v), causes no change in the 
velocity profile. These considerations can easily bring down the 
true values of P^ for our calculations well below the approxi- 
mations quoted above. In any case, there is currently no general 
consensus about the correct value for P^: we note that rather 
high values have recently been found in some numerical MRI 
simulations (e.g. see Takahashi & Masada 120101 Romanova et 
aI. l20TlT l. 
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(a) 770 = 4 X 10'" cm^ s"' 
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Fig. 7: Poloidal magnetic field lines. Comparison between two Fig. 8: Contour plots of B^. Comparison between two configura- 
configurations having values of 770 larger than that of the refer- tions having values of 770 larger than that of the reference config- 
ence configuration. uration. 



We chose our values for 77 bearing in mind which values 
would give significant field distortions in the disc and there- 
fore be most interesting. However, it is clearly important to in- 
vestigate the effects of varying these numbers, and so we also 
made some calculations using larger values of 77 (smaller P^). 
We show here results for 770 = 4 x 10'" cm^ s"' and 770 - 10" 
cm^ s"' with 771; being, as usual, two orders of magnitude larger 
(in the reference configuration we used 770 = 10"^ cm^ s"'). 

Results for both configurations are presented in Figs. |7] [8] 
and|9l where we show the poloidal field lines, the contour plot 
of the toroidal field component and the radial profile of its shell 
average, respectiveljQ. 

Increasing tjo (i.e. decreasing the magnetic distortion func- 
tion Dm, our generalisation of the magnetic Reynolds number 
introduced in Paper I) makes the field diffuse more efficiently 
and therefore the solution gets progressively nearer to being a 
dipole field (see Fig. |7]l. Reducing the poloidal field distortions 
in turn changes the toroidal component, and the modifications 
are quite substantial. The structure with four extrema gradually 
turns into one with only two (see Fig.jSJ, where the sign of is 
positive for radii smaller than the corotation radius and negative 
for larger ones. This transition is clearly seen when plotting the 
shell average of the toroidal field (see Fig.|9]l. The behaviour of 
the sign of the shell average of is then the same as that for B^ 



^ We have also run a configuration with 770 = 10'^ cm^ s ' and found 
results completely in agreement with the trends shown here. 
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Fig. 9: Shell average of the toroidal component. Comparison be- 
tween two configurations having turbulent magnetic diffusivity, 
770, larger than that of the reference configuration. The values of 
770 are 4 X 10'" cm^ s"' for the solid curve and 10" cm^ s"' for 
the dashed curve. 



in the early analytic models, although the details of the profiles 
have significant differences (see Fig.[T3]below). 

Finally we note that the magnitude of the toroidal component 
decreases with increasing 770, the maximum of the shell average 
being ~ 14x lO"-' Bo for the configuration with 770 = 4x lO'" cm^ 
s"' and ~ 6 X 10"^ Bo for the one with 770 = 10" cm^ s"' (com- 
pare the shell averages in Fig.|9]with each other and with the one 
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Fig. 11: Shell average of the toroidal component for 77 = 10" 
cm^ s"'. For // - 10^^ cm^ s"', the curve has exactly the same 
shape, but the values are about 100 times smaller (in absolute 
value). 



in Fig. |6]l, while the maximum for the poloidal component has 
remained at ~ 3.5 x lO^^* Bo (this is because the maximum for 
the poloidal componant occurs at the inner edge, where the field 
depends more on the boundary conditions than on the value of 
ri). 



4.3. Configurations witli constant ij 

Although the rj profile that we have used so far is the one that 
we consider to be the most appropriate when studying accre- 
tion within the kinematic approximation (for the reasons given 
in Sect. |2]i, we have considered also configurations where the 
difFusivity is constant through all of the disc and the corona, in 
order to have a clear understanding of how sensitive the results 
are to this quantity. Cases with 77 = 10'" cm^ s"' and 4x 10'" cm^ 
s"' proved unsatisfactory because of having very abrupt changes 
away from the dipole field at the top of the corona (this is exactly 
the behaviour that we have tried to avoid by using a larger value 
of T] in the corona in our reference configuration). Results for 
77 - 10" cm^ s"' and t] - 10'^ cm^ s"' are shown in Figs. fTOlfTTI 
andO 

For configurations with constant 77, the analysis is made sim- 
pler since there are no regions with diffusivity gradients. As 
stated in the previous subsection, using a larger value of 77 re- 
duces the deviations away from the dipolar field. However, in 
contrast with the previous case, the distortions are now more uni- 
form throughout the disc (compare Figs. |7] and [TOli. because 77 is 
constant and the magnetic distortion function is monotonically 
decreasing with r (while previously drD^ had a peak). We note 
that deviations away from a pure dipole are very small when us- 
ing 77 = 10'^ cm^ s"' , and for 77 comparable to or larger than this, 
the poloidal component of the field can be well approximated 
by a dipole. As regards the toroidal component, it has only two 
extrema which are both located just beneath the surface of the 
disc (see Figs. [TT] and [T2l i. These two extrema have been ob- 
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served in all of the configurations which we have studied; they 
are also present in the models of Wang (1987) and Campbell 
( fT987.,1992j and therefore seem to be robust features. 

5. Discussion 

5.1. Comparison with analytic modeis 

In the analytic models of Wang (1 19871 1 and Campbell (I1987I I. 
who we will refer to now as W&C, the toroidal component of 
the magnetic field is written as being proportional to the relative 
angular velocity between the disc and the star multiplied by the 
vertical field, which is taken to be a pure dipole (see Eq. ([TJ). For 
models where Qdisc is Keplerian and the inner edge of the disc 
is inward of the corotation point, has a positive global maxi- 
mum at the inner edge of the disc, is zero at the corotation point, 
reaches a global negative minimum somewhere outward of this 
and then tends towards zero at very large r. If Q. deviates from 
Keplerian in the inner part of the disc (reaching a maximum and 
then decreasing again as one moves inwards), then the location 
of the maximum of is not in general at the inner edge, but 
depends on the precise profile of O in this inner region. 

The above description is only partially in agreement with the 
results of our present two dimensional calculations. They share 
the feature of having a zero of B^ at the corotation point (or ex- 
tremely close to it, see Fig. |6]). As regards the predicted global 
maximum of B^ in the inner part of the disc, all of our calcu- 
lations show a positive maximum close to where Q. has a maxi- 
mum (compare Figs. l5l [8] and [T2l with Fig. |2]l. However, when 
77 is not constant, magnetic field lines accumulate and a sec- 
ond maximum can appear outward of the corotation point, and 
this can also become a (positive) global maximum depending 
on the value of 770. Finally, regarding the minimum: as in the 
W&C models, we always find a global negative minimum be- 
fore the outer edge of the disc; however, when rj is not constant 
a second minimum can appear, thus producing a structure with 
four extrema: two maxima and two minima (see Fig.Hl). We note 
here that even in the configurations with only two extrema, as in 
W&C, the profile of the toroidal component is still different from 
that predicted by those models. In particular, the location of the 
minimum and the magnitude of the field at both extrema can be 
very different. 

The main features of the toroidal field component, as given 
by our present numerical calculations and by the analytic mod- 
els, can be seen in Fig.[T3j where we show three profiles for B^: 
the dotted curve is the W&C profile for our model, as resulting 
from Eq. ([T]i but where we have used our profile for the angular 
velocity near to the equatorial plane (i.e. Eq. ( |23] ) at = 88°), the 
solid curve shows the shell average for our reference configura- 
tion, and the dashed curve is for the configuration with constant 
77(101' cm^ s-i). 

The properties of having additional extrema of B^ and of 
varying the location and magnitude of the two standard extrema, 
are not seen in the W&C models. 

5.2. Role of Bp and 

In Fig. [14] we show the quantity AQ Bg, where AO = Odisc - f^star 
and Bg is the 0-component of the poloidal field as obtained from 
our numerical calculations for the reference configuration. This 
quantity has one global maximum and one global minimum. The 
maximum is located very close to where B^ and Q have their first 
maxima, and the minimum is at a radial location close to that of 
the first minimum of B^ (only a few smaller - see Fig. |5] and 
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Fig. 13: Comparison of three different profiles of the toroidal 
field: (1) the solid curve is the shell average for the reference 
configuration, (2) the dashed curve is the equivalent one for the 
configuration with constant 77 (10" cm^ s"'), (3) the dotted curve 
is given by the W&C formula with our Q profile at = 88°. All 
of the curves have been normalised so as to have their peak at 1 . 
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Fig. 14: Contour plot of AQ Bg for the reference configuration. 



TableO. This shows that the quantity AQ.Bg is still relevant for 
grasping the fundamental properties of the toroidal component 
of the field, although care must be taken in choosing the profile 
of Bg. Some differences between the predictions of the W&C 
models and our numerical results can, in fact, be explained in 
terms of the distortion of the 6 component of the magnetic field. 
However, in Fig. [14] there is no evidence for the additional max- 
imum and minimum, and so this is not the whole story. We need 
to consider the distortions of the field in more detail, and not 
focus only on the quantity AQ Bg. 

In Paper I we have shown that the distortions of the poloidal 
component due to the plasma motion can be well described by a 
generalisation of the magnetic Reynolds number (which we have 
called the magnetic distortion function), defined afl 

Du. = — . (30) 



Note that we have included vg here in the definition of (it was 
not present in Paper I, because we were showing there results on the 
equatorial plane). We recall that vq is zero in the disc and equal to 
0.5 V,. tan Q in the corona. 
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Fig. 15: Radial derivative of the magnetic distortion function 
near the equatorial plane (at 9 - 89°) for three configurations: 
(1) the solid line is for the reference configuration with rjo - 
10^° cm^ s"', (2) the dotted line is for the configuration with 

7;o = 410l°' 



' cm^ s 



, (3) the dashed line is for tjq = 10" cm^ s 



We studied Dm in detail in Paper 1 (see section 4 of that pa- 
per); here we just recall that the magnitude of the peak in its 
radial derivative is a measure of the degree of accumulation of 
poloidal field lines in its vicinity. For the cases with constant 
the field amplification caused by this distortion is negligible 
(the distortion is more homogeneous and the field lines do not 
accumulate) and is just proportional to v, in all of the disc 
(so that there is no peak at all in the radial derivative). For the 
reference case, instead, 5, does have a peak and its radial lo- 
cation (r ~ 28.5 Tg) is near to that of the additional maximum 
(r ~ 30 Tg). Increasing tjq reduces the magnitude of the peak and 
also the additional maximum of becomes less pronounced, 
eventually disappearing for 5, Dm < 0.037 (see Fig.fTSl). We can 
therefore draw the conclusion that it is the radial derivative of 
the magnetic distortion function which is responsible for the ad- 
ditional maximum in the toroidal field profile. 

5.3. Our picture for 

In this subsection, we develop our alternative picture for the 
structure of the toroidal component of the magnetic field, fol- 
lowing the same general approach as W&C, but with a more 
detailed representation of the poloidal magnetic field and the ve- 
locity field, and a two-dimensional model. This is still a very 
simplified picture but we believe that it can represent a useful 
step forward, giving some additional insights. Our starting point 
is the (p component of the induction equation: 



drB^ = [V X (V X B)]^ - [V X (77V X B)]^ 



(31) 



The advective term is the one responsible for generating the 
field, while the diff'usive term causes field losses. We then de- 
fine the following scalar quantities: 



G = d,B^W = [Vx(vxB)]^ 
L = d,B^U = [V X (77V X B)], 



(32) 
(33) 



In a steady state dtB^\+ - d,B^\- so that G - L and d,B^ = 0. 

We focus first on the generation term and split it into two 
parts: one involving the poloidal motion (Gp) and the other in- 
volving the toroidal motion (G^): 




Fig. 16: Contourplot of ao in Cartesian coordinates, in the region 
r = [10 - 115] rg. This should be compared with Fig.H] 

In spherical coordinates, these two terms are written as 

Gp = [ V X (Vpo, X B)]^ = - i [drir V, B^) + deivg B^)] (35) 

G^ = [Vx (v^ X B)]^ = i [dAr B,) + dg{v^ Bg)] . (36) 

From these expressions we can see that, in general, the genera- 
tion rate depends on several quantities and not only on the ver- 
tical gradient of the angular velocity (as in the W&C models). 
In fact, all of the components of the velocity field and magnetic 
field are present. However, usually in discs » vv, vj^and one 
can expect the second term to dominate. An important diff'er- 
ence between Gp and G^ is that the former cannot generate any 
toroidal field from zero, but can only modify B^ once it has al- 
ready been produced by some other means, e.g. coming from 
advection by the azimuthal motion (i.e. from G^). 

We have calculated the ratio G^/Gp for our configurations 
and have found that the toroidal term always dominates, even 
if the ratio is not as large as would have been expected just by 
comparing the velocities (because one should consider also the 
magnetic field). However even neglecting the contribution from 
the poloidal advection, the generation rate for B^ is still consid- 
erably different from the one considered in the W&C models. In 
fact, it contains also the radial derivative of the term r vv Br and 
the vertical derivative of Bg. 

So far, we have focused only on the generation term. In order 
to obtain the profile of the toroidal field we have to equate it to 
the loss term, for which we have 



L = [V X (77V X B)]^ 

= --U[T]dr{rB^)]+de 



1 



r sin0 



dg [b^ sin 6*) 



(37) 
(38) 



If we consider only typical values {rj and r), then we can write 

(39) 



In a steady state, when L - G - G^, one has G^ - (fi/f^)B^, and 
one gets 



Bd, 



Gib — . 
1 



(40) 



dtB^\+ — Gp + Gfj, 



(34) 



' In some circumstances one could have a strong wind from the top 
surface of the corona and v,. « vg < v^. 
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Fig. 17: Contour plot of Q for the reference configuration. All 
locations with Q > I are marked in red so as to highlight the 
regions where Q < 1. 



In our configuration, has two main characteristic values: rjo in 
the main disc region and ?/c - 10' rjQ in the corona and in the 
inner part of the disc. If in Eq. ( l40l l we replace rj with the actual 
profile of the magnetic diffusivity, and choose r - rg, then we 
find that the variation of Bj, can be well approximated by that of 
the coefficient ao in Eq. ([8J: 



rrj 



dr(r V0 Br) + deiy^ Bg) 



Bq ao 



(41) 



where Bq is the reference unit of the magnetic field and we recall 
that in Eq. ^ the fields are dimensionless (while here they are 
dimensional). 

This formula for B^ is rather well confirmed by our numeri- 
cal calculations. In Fig. [16] we show the contour plot for a^: this 
is very similar to that for the toroidal field (compare with Fig.lUi, 
the differences being due to the approximations made in evaluat- 
ing the loss term L, and to neglecting Gp in the generation term. 

5.4. Consistency check with the W&C models 

From Eq. (1411 1 we can see that the profile of the toroidal field is 
basically determined by four factors: (1) the radial derivative of 
Br, (2) the radial derivative of rv^, (3) the vertical derivative of 
Bg and (4) the vertical derivative of v^. 

In the W&C models, the first three of these are neglected 
because Br is taken to be zero everywhere and Bg is supposed 
not to vary with 0. If one adds also the other assumptions used in 
their models (about the velocity profile and the disc thickness), 
one then finds that Eq. (HTt reduces to 



B4, — {^K - ^s) Bg — 

n 77 



(42) 



which is the same as Eq. ([TJ with -/„ - ("nr/h) and - (rg/j]). 
We can therefore recover the expression appearing in the earlier 
analytic models from our result. 

We now want to check a posteriori whether or not the sim- 
plifications made in the W&C models are still valid in our more 
general 2D model and, if so, in which parts of the disc. One can 
do this by calculating the ratio (Q) between the two terms on the 
right-hand- side of Eq. dTO : 



\dg(v^Bg)\ 
\dr(rv^Br)\ 



Clearly, a necessary condition required for matching with the 
W&C models is that 2 » 1. In Fig.fTTlwe show a contourplot of 
this quantity for the reference configuration. We show in red all 
of the regions where Q > I, while regions having the condition 
clearly being violated are colour-coded to show by how much Q 
is smaller than 1 . 

There are large parts of the disc where this necessary con- 
dition is not met, with Q even reaching values as small as 10"^. 
Moreover, even in regions where Q » I, the further necessary 
condition dgBg <K dgO. may not be met. In the main disc region 
the angular velocity is almost constant with 6 (the transition from 
Keplerian rotation to corotation happens in the corona), and here 
we are exactly in an opposite regime, i.e. dgBg » dgQ. ~ 0. Only 
in the corona and in the upper part of the disc are the vertical 
gradients of Q not negligible. 

For a pure dipolar field, Eq. (l43T l becomes 



Q 



dip 



I-Hicote^ 



-1 



(44) 



which, for Keplerian motion, gives exactly 2/5. Therefore the 
necessary condition is never met for a pure dipole and Keplerian 
rotation, which is a good description for the parts of our discs 
near to the mid plane. This should not surprise us, since two 
key assumptions made in the W&C models hold only in differ- 
ent parts of our discs and not together: the field was supposed 
to have both B, and dgBg vanishing and decaying as for a dipole 
(which we have only very close to the equatorial plane or for 
large //) and the transition to corotation in the angular veloc- 
ity was supposed to be very sharp (which we have just beneath 
the disc surface, where Q. has to become equal to Qj, far from 
the equatorial plane). We have calculated Q for a configuration 
with constant t] (10'^ cm^ s"'): in this case the diffusivity is so 
strong that deviations away from the dipole are quite small and 
so changes in Q come almost entirely from the angular velocity 
(according to Eq. (l44l i). As expected, in the lower part of the disc 
we find Q ~ 0.4. 

Therefore having a larger value of 77 does not necessarily im- 
ply that the necessary conditions hold in a larger region of the 
disc, it just implies that the field is closer to a dipole (which is 
what we are imposing at the boundaries). In order to get Q to go 
to infinity (which is what was assumed in the W&C models) not 
only does one need Br -0 and dgBg - 00 but also that the verti- 
cal gradient of the angular velocity, at the same location, has to 
be non-zero, or much larger than the vertical derivative of Bg. 

6. Conclusions 

In this paper we have considered a system consisting of a ro- 
tating neutron star, having a dipole magnetic field aligned with 
the rotation axis and surrounded by an accretion disc. The disc 
is truncated at the Alfven radius and has a coronal layer above 
and below it. The region outside the corona is taken to be vac- 
uum and we impose dipolar boundary conditions on all of the 
boundaries. 

Our aim was to improve on the analytic models developed 
by Wang ( 1987 ) and Campbell (fT987] l (W&C). As in those mod- 
els, we have made the kinematic approximation and have looked 
for an axisymmetric stationary solution of the induction equa- 
tion, but we have gone beyond those models in solving for all 
of the components of the magnetic field and not assuming the 



(43) 



* These conditions hold exactly for a dipolar field only at the equato- 
rial plane; as one moves away from that, they are only partially satisfied. 
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poloidal component to be dipolar within the disc. We have also 
retained all of the components of the velocity field rather than 
putting Vr and vg to zero everywhere. We have performed a fully 
two-dimensional calculation, without making any vertical aver- 
age or Taylor expansion in h (the semi-height of the disc). Finally 
we have neglected dynamo action but have included a turbulent 
magnetic diffusivity. 

The analysis of the poloidal component of the magnetic field 
has been presented in a previous paper (Naso & Miller, 120101 
Paper I); in the present paper we have focused on the toroidal 
component. We have solved the (p component of the induction 
equation numerically and have shown that the profile obtained 
for can be very different from that in the earlier analytic mod- 
els. 

In the W&C models, the toroidal field strength was taken to 
be proportional to the relative angular velocity between the disc 
and the central object multiplied by the vertical field, which was 
taken to be dipolar. However in Paper 1 we found that, when cal- 
culated consistently, the poloidal field component was often far 
from being dipolar. Therefore a first improvement with respect 
to the earlier models was to use the poloidal field as obtained 
in our calculations, i.e. a field dragged inwards by the plasma 
motion. This behaviour explains why we then find different in- 
tensities for the toroidal field, and also a different location for its 
global minimum. 

When the turbulent magnetic diffusivity rf increases or the ra- 
dial velocity decreases one expects the field to be progressively 
less distorted by the plasma motion. This is indeed what we have 
found both here and in Paper 1. Our results show that when the 
diffusivity 770 is larger than about 10'^ cm^ s"', with the charac- 
teristic velocity vq being of the order of 10^ cm s"' , then the field 
is barely modified. Therefore whenever we expect the magnetic 
field to deviate from the stellar dipole, 77 should not be larger than 
10^ Ivcgsl cm^ s"', where |vcg,sl is the characteristic magnitude of 
the radial velocity expressed in cm s"' . 

When the turbulent magnetic diffusivity is not constant 
throughout the disc and corona, two additional extremal points 
may well appear: if the radial derivative of the magnetic distor- 
tion function is larger than a critical value (about 0.04 in the 
equatorial plane), there is an additional maximum and minimum, 
and in some cases can even become positive again outward of 
the corotation point, so that there are additional locations where 

= 0. It is clear that under these circumstances the picture of 
which regions of the disc tend to spin the star up or down has 
to be radically redrawn (this will be the subject of a future in- 
vestigation). However we should emphasise here that there are 
still many uncertainties among experts about which profile of rf 
should be used and we have therefore made very simple choices 
here in line with our step-by-step approach. 

We have presented a new suggestion for the profile, which 
reduces to that of W&C if one imposes B,- - dgBg - and 
Q = Q.K- In general there are large parts of the disc where the ad- 
ditional terms included in our new picture for B^ dominate over 
the one retained by W&C (see Fig. [TTl i. Our simplified expres- 
sion (Eq. |4TI ) reproduces the numerical results quite well (com- 
pare Figs. l4l and [TSIl, the differences being due to approximations 
made in calculating the generation and loss terms for B^. 

Summarising, in our calculations we have found that B^ can 
have two maxima and two minima (see Fig. |4li. The first maxi- 
mum (positive and inward of the corotation point) and the first 
minimum (negative and outward of the corotation point) can be 
explained referring to the quantity AfiBe, which has two ex- 
trema at the same locations as for B^ (see Fig. [T4l i. These ex- 
trema appear also in the W&C models, where the toroidal field 



is, in fact, taken to be proportional to AQ B,. There is a funda- 
mental difference however: in W&C B- is a pure dipole, whereas 
the Bg which we consider here is that of a field being dragged in- 
ward by the motion of the accreting material. When rj is not con- 
stant, there is an additional maximum whose magnitude and sign 
depend on the diffusivity in the disc, and whose radial location 
is always outward of the first minimum, coinciding with that of 
the maximum in the radial derivative of the magnetic distortion 
function D^. Outward of this maximum, the field tends to come 
back to the profile that it would have had if 77 were constant, and 
this produces the last minimum (compare Figs. l9l and [TTI). 

The main conclusion of this analysis is that, when the 
poloidal component of the magnetic field is treated self- 
consistently in the calculations, the profile for B^ can be sig- 
nificantly different from that obtained by W&C, and the mag- 
netic torque generated by it would then be different as well. 
Moreover, when the turbulent magnetic diffusivity rj is not con- 
stant throughout the disc and corona, some additional unex- 
pected features can appear (such as a region of positive B^ out- 
ward of the corotation point). In the present work, we have re- 
tained the very simple Keplerian rotation law in the main part of 
the disc. Even within a purely hydrodynamical treatment, more 
complicated velocity fields than this are expected (see Kluzniak 
and Kita, 120001 Jiao and Wu, 1201 II ) and further changes are ex- 
pected when back-reaction from the magnetic field on the veloc- 
ity field is included. The effects of this will be another topic for 
investigation in subsequent stages of our step-by-step approach. 

Appendix A: Testing of the code 

In this Appendix we discuss some of the tests that we have per- 
formed on the numerical code used to solve Eq. ([8]l. For a de- 
scription of the Gauss-Seidel relaxation procedure and of the 
discretization scheme see Appendix A. 1 of Paper I. 

Before describing the tests, we should underline a differ- 
ence in the boundary conditions with respect to the code used 
for the poloidal analysis. In Paper 1 we were not imposing the 
dipolar boundary conditions on the magnetic field directly but 
rather on the magnetic stream function S. The poloidal mag- 
netic field was then calculated by differentiating S. Because of 
this, B, and Bg were not precisely dipolar on the boundaries. For 
the calculations in the present paper, we have introduced a row of 
ghost points, where we set the field to be exactly dipolar, i.e. the 
poloidal component is a pure dipole and the toroidal component 
is zero. 

We divide the tests into two groups. In the first group we 
chose the coefficients of Eq. ([8]) in such a way that it was possible 
to find an analytic solution, while in the second group we used 
the values given by Eqs. (l9])-(fT3b. Here is a schematic description 
of these tests: 

1 . Tests with analytic solutions 

In addition to fixing the coefficients, one also has to choose 
the boundary conditions. We considered three sub-cases: 

1 . 1 All coefficients constant and set to 1 .0. There is then the 
following analytic solution: 

B^ = exp(-x/10^ - 012) cos ( V2.96046I/2) - 1 . (A.l) 

1.2 All coefficients set to 0.0 except for 

a„ = 1 Ix^ (A.2) 

a^^blx. (A.3) 
The general solution is then: 

B^ = /i exp (/ 0) ( A.4) 
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where k\ is a function of b and k2. We chose h - 10, 
k2 = 36^ and b - 2 - kj, so as to include a complete pe- 
riod of the angular part of the solution within our angular 
domain. The solution is then: 

B^^ — cos (36 9) . (A.5) 

X 

1.3 The same choice of coefficients as in test 1.2 but with 
different boundary conditions: h - \Q, k2 - Q and b - 
70. The analytic solution is then: 

= 10 . (A.6) 

2. Testing the model setup 

In this group of tests we used a setup which was very similar 
to that used for our actual physical analysis but varied some 
numerical parameters so as to test the code. We performed 
four tests aimed at: 

2.1 checking convergence by changing the number of grid- 
points; 

2.2 studying the dependence of the solution on the location 
of the radial outer boundary; 

2.3 studying dependence on the initial estimate for the solu- 
tion; 

2.4 optimising the iteration step size by simple benchmark- 
ing. 

A. 1 . Tests with analytic solutions 

For all of the tests in this category, we considered the code sta- 
bility and convergence. We used grids with different numbers of 
points in both directions and compared the analytic errors and 
the solutions. 

In all cases, the stability of the code was related to the size of 
the iteration step, the code being stable for values smaller than a 
certain threshold. 

To check convergence, we considered how the maximum of 
the analytic errors changed with varying the total number of iter- 
ations. We considered both absolute and relative errors and anal- 
ysed them in the region of interest (i.e. for r < ric) by calculating 
their maximum and looking at their 2D profile. As the iteration 
procedure continues, the errors decrease and at a certain point 
they saturate, so that making more iterations no longer leads to 
smaller errors. In addition to the errors, we have also considered 
the evolution of the root mean square (rms) of the solution. 

The size of the errors at saturation depends on the grid reso- 
lution, being smaller for grids with more points (the size of the 
domain is fixed, so that increasing the number of points means 
increasing the resolution). Moreover the improvement obtained 
when increasing the number of grid-points becomes progres- 
sively smaller, as it should do in a convergent regime. We cal- 
culated an effective order of convergence peft by considering 
the maximum value of the relative error at the final iteration, 
and how it changed with the grid size. Doing this we obtained 
Peff -1.5. 

All of the tests gave satisfactory results, confirming stabil- 
ity and giving convergence within ~ 10"* iterations for tests 1.1 
and 1.2, and within ~ 10^ iterations for test 1.3. The maximum 
saturation error was ~ 10"^% for test 1.3 and ~ 10"^% for test 
1.1. As regards test 1.2, since the solution has some zeros in the 
considered domain, we could not estimate the error by consider- 
ing the relative error (since it diverges at the zeros). We instead 
considered the maximum of the absolute error and compared it 
with the rms of the solution, finding that it was less than 1%. 



A.2. Tests with the model setup 

For these tests we had no analytic solutions, and so could not 
calculate analytic errors but only residuals. When considering 
stability, we looked at the residuals, while for testing conver- 
gence we considered the change in the 0-average of the solution 
(or in its rms) when changing the number of grid-points. 

We recall that we cannot change the number of grid-points 
freely. In fact, one of the coefficients (aq) is not calculated from 
an analytic function but comes instead from the numerical re- 
sults of Paper I for the poloidal field. This coefficient is therefore 
directly defined only on the grid used in that calculation, which 
had 1001 X 21 points. We consider this grid as our reference one. 
When using a grid with fewer points, we have to choose them as 
a subset of our reference grid, while when using a larger number 
of points, we have to interpolate. Similarly when changing the 
size of the domain (by reducing rout), we must also decrease the 
number of grid points A^,- accordingly, as explained at the end of 
Sect.llKsee also Eq. (E9]l). 

A very important result of test 2.1 is that if the transition in 
77 between the disc and the corona is not well-enough resolved, 
then we find convergence to a different solution. More precisely, 
using an angular width for the transition of 5 x 10"^ radians we 
need to have at least 20 grid points overall in the 6 direction in 
order to converge to the correct solution, i.e. to the same solution 
as for grids with larger values of A^^ (we tested with A^^ = 39 
and Nj - 77). This gives a minimum number of angular zones 
required for convergence of about 5 within the transition region. 

Test 2.2 tells us that the solution is not very dependent on the 
location of the radial outer boundary, contrary to the situation in 
Paper I for the poloidal field. The ^-averaged solution obtained 
using rout = 290 differs from that obtained with rout = 750 by 
less than 0.02% in all of the region of interest. If we consider 
the solution rms, the difference is even smaller, being less than 

5 X 10-'*%. 

For test 2.3, we ran calculations with quite different, and 
even unphysical, initial estimates for B^, in order to see whether 
the solution still converged correctly. We used (1) a constant 
value B^ - 1; two decaying profiles: (2) B^ - (10/x)^ and 
(3) B^ = (10/x); (4) a growing profile B^ = Vjc/TO and (5) 
a Gaussian profile in both directions (centred at x = 100 r^, 

6 = 85° with widths 15 r^ and 2°). We obtained the same fi- 
nal solution for all of the cases; profiles (3) and (5) converged 
after ~ 2 x 10^ iterations, while all of the others converged after 
~ 4 X 10* iterations. 

Finally in test 2.4 we changed the iteration step size Af, look- 
ing for the largest possible value still giving stability. As in Paper 
I, we found that the maximum value depends more sensitively on 
Nj than on A^,. For A^j = 81 we found Aw = 4.7476 x 10""*. 
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